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ABSTRACT 

We consider the thermal structure and radii of strongly irradiated gas giant planets over a range in 
mass and irradiating flux. The cooling rate of the planet is sensitive to the surface boundary condition, 
which depends on the detailed manner in which starlight is absorbed and energy redistributed by fluid 
motion. We parametrize these effects by imposing an isothermal boundary condition T = Tdeep below 
the photosphere, and then constrain Td ccp from the observed masses and radii. We compute the 
dependence of luminosity and core temperature on mass, Td ee p and core entropy, finding that simple 
scalings apply over most of the relevant parameter space. These scalings yield analytic cooling models 
which exhibit power-law behavior in the observable age range 0.1 — 10 Gyr, and are confirmed by time- 
dependent cooling calculations. We compare our model to the radii of observed transiting planets, 
and derive constraints on Tdeep- Only HD 209458 has a sufficiently accurate radius measurement 
that Tdeep is tightly constrained; the lower error bar on the radii for other planets is consistent with 
no irradiation. More accurate radius and age measurements will allow for a determination of the 
correlation of Td ee p with the equilibrium temperature, informing us about both the greenhouse effect 
and day-night asymmetries. 

Subject headings: planetary systems-planets and satellites:general 



1. INTRODUCTION 

Foll owing the discovery of the planet orbit ing 51 
Peg lIMavor fc Ouelo^l QflStlMaTc'v fc Butlerll 995l). more 
than 160 planets have been found around nearby stars us- 
ing precision Doppler spectroscopy. 1 Theories of planet 
formation now have the demanding task of explaining the 
existence of gas giants with semi-major axes one hundred 
times smaller than Jupiter, others with order unity or- 
bital eccentricities, a detailed spectrum of (minimum) 
planet masses, and metallicity correlations with the par- 
ent star. 

The discovery of transiting planets in the last five years 
(see Table nj challenges not only theories for the origin 
of short-period gas giants, but also their structure and 
thermal evolution, spectrum, and interior fluid dynamics. 
Measurements of planetary mass, radius, and (stellar) 
age test cooling models which predict radius as a func- 
tion of mass and age. The atmospheres of two planets 
have been directly observed. For HD 209458b, absorp- 
tion lines (due to stellar photo ns passing through planet's 
atmosphere) have been fo u nd itCharbonneau et alJl2002t 
iVidal-Madiar et all 120031 I2004[K and the first detec- 
tions of photons emitted by planets outside our so- 
lar system have bee n made for the the rmal emission 
from HD 209458b jDeming et alJ 12005^) and TrES-1 
ijCharbonneau et alJl2005|) . These observations directly 
constrain the atmospheric structure, temperature profile 
and chemical composition near the photosphere. 

Evolution of the short orbital period transiting exo- 
planets is significantly different than for Ju piter and Sat- 
urn d ue to proximity of the parent star (Guillot et al. 
1996). Irradiation increases the photospheric temper- 
ature by nearly an order of magnitude relative to an 

1 For up to date catalogs, see [^ ^g^//e^ojgjaiie ts.org/| and 
http: / /obswww.unige.ch~udry/planet/planet.html 



isolated planet. Irradiation also decreases the cooling 
rate, and hence the rate of shrinkage, by alteri ng the 
surface boundary condition (B urrows et al.lf2000|) . This 
is immediately apparent in Table ^ as many transiting 
extra-solar giant planets (EGP's) have radii significantly 
larger than Jupiter. As short per iod planets are ex- 
pected to be tida lly synchronized l)Ouillot et alJ Il996t 
iMarcv et alJ Il997). the strong day- night temperature 
contrast will driv e winds to transport heat f rom the day 
to the night side ({Showman fc Guillotll2002T) . Hence the 
atmospheric temperature profile depends on a combina- 
tion of detailed radiative transfer calculations for absorp- 
tion of starlight, and hydrodynamics to model day- night 
winds and dissipation of wind kinetic energy. Lastly, 
tides raised on the planet by the par ent star may signif- 
icantl y affect its thermal evolution ijBodenheimer et alJ 
l2pj0lt IShowman fc Guill<^l2W)2]K The free energy avail- 
able by synchronizing the planet's spin or circularizing 
the orbit are comparable or larger than the thermal en- 
ergy. Hence if the heat can be deposited sufficiently deep 
in the planet in less than a cooling timescale, the cooling 
can be slowed, or even reversed. However, i t is uncertain 
if tides can deposit heat deep in the planet l|Lubow et alJ 
lL997tlOgilvie fc Linll200l IWJI2004D . 

Evolutionary models show that the cooling rate is 
quite s ensitive to the uncertain surface boundary con- 
dition ijGuillot fc Showman! 12002]) . This boundary con- 
dition has been implemented using vario us approxima- 
tions . Full radiative trans fer calculations l)Barman et alJ 
120011: iHubenv et al.ll2003|) of static atmospheres include 
the stellar irradiation self-consistently, and determine 
the temperature structure for a given cooling flux from 
the deep interior. These calculations compute (rather 
than assume) the albedo, and determine the tempera- 
ture rise due to absorption of starlight (the greenhouse 
effect). Such detailed radiative transfer solutions have 
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TABLE 1 
Transiting Extrasolar Planets 



object o(au) M p (Mj) R P (Rj) T C<1 [K] a T dccp [K] h Age (Gyr) Reference 



OGLE-TR-132 


0.031 


1.19 ±0.13 


1.13 ±0.08 


2100 


< 2200 


0-1.4 


1 


OGLE-TR-56 


0.023 


1.24 ±0.13 


1.25 ±0.08 


2100 


1000 - 3100 


3±1 


2,3,12 


HD 209458 


0.046 


0.69 ±0.05 


1 ni+0.05 

L - 3L -0.05 


1500 


2200-2800 


4-7 


4,5 


OGLE-TR-10 


0.042 


0.63 ±0.14 


1.14 ±0.09 


1500 


< 2600 




6,12 


OGLE-TR-113 


0.023 


1.35 ±0.22 


1 n8+ - 07 


1300 


< 2100 




7 


TrES-1 


0.039 


0.73 ± 0.04 


1 ns+ - 05 


1200 


< 1000 


2.5 ± 1.5 


5,8 


OGLE-TR-111 


0.047 


0.52 ±0.13 


0.97 ±0.06 


1000 


< 1200 




9,12 


HD 149026 c 


0.042 


0.36 ± 0.04 


0.725 ±0.05 


1700 




2.0 ±0.8 


10 


HD 189733 


0.031 


1.15 ±0.04 


1.26 ±0.03 


1200 


< 3200 




11 



References. — (1) IMoutou et~aH ll200l) (2) iTorres et alJ l20bH) . (3) Sassolov (2003) , (4) 
ICodv fc Sasselovl (120021). C5) iLaughlin et alJ 120051). C6 ) lKonacki et alJ (1200,^ . f7)| Bouchv et al. 
l|2004ft . f8) ISozzetti et al1ll2004l) . f9) lPont et al.ll|2004D . (TO) ISato et alJ lj200fil ) .(ll lBouchv et al. 
pOOa) . Q2) ISantos et all 12000) 

a Here T oq = T t (R t /2a) 1/2 . See the discussion following eq. Q. 

b Allowed range of Td ccp given range of mass, radius and age. If no age range given in the 
literature, we (arbitrarily) give the maximum value of Td eep for an age less than lOGyr. However, 
given an accurate age range, the figures in §|H]can be used to obtain stronger constraints than 
given here. 

C HD 149026's small radius clearly indicates a large core size or heavy element abundance. The 
present paper does not include heavy element cores, so we do not discuss HD 149026 further. 



been incorporated as boundary conditions for some evo- 
lution ary calculations l)Baraffe et al!2f)03tlBurrows et all 
2003). However, as day-night and equator-pole winds are 
not included, assumptions must be made about how the 
stellar flux is deposited over the surface of the planet 
(only day-side versus evenly over the entire surface, etc.) 
which directly affect the temperature profile. Other evo- 
lutionary calculations (e.g. Bodenheimer et al. 2003) 
solve the radiation diffusion equation and set the tem- 
perature at (infrared) optical depth 2/3 to be the equi- 
librium temperature, ignoring additional temperature in- 
crease due to absorption of starlight. Lastly, a num- 
ber of groups llShowman fc Guillotj2002HCho et al.l2003t 
[ Burkert et alll2005l iCooper fc Showmanll2005l llro et all 
2005) are beginning to model the day-night winds on 
tidally lo cked, short orb i tal pe riod planets, and the role 
of clouds iFortnev et all l)2003|) . As we stress here, the 
crucial parameter for the cooling rate is the temperature 
at the radiative-convective boundary, which is orders of 
magnitude deeper in pressure than the photosphere. 

The plan of the paper is as follows. The uncertain sur- 
face boundary condition is discussed in § El motivating 
the surface isotherm used in our models. Details of cool- 
ing models and microphysical input are described in § 
In § 0] we compute the dependence of the luminosity on 
planet mass, core entropy and irradiation. An analytic 
solution for the temperature profile in the radiative zone, 
and the position of the radiative-convective boundary are 
derived in § These results are collected together in § 
to derive an analytic cooling model which exhibits simple 
power-law dependence on time. The radii of irradiated 
gas giant planets are discussed in § and the analytic 
formula for the radius given in eq. (|32[) . We apply our 
models to the observed transiting planets and give con- 
straints on the temperature of the deep surface isotherm 
in § |H1 Our main conclusions are summarized in § El 



2. SURFACE BOUNDARY CONDITION 

The surface boundary condition we adopt is to set the 
temperature T = Tdeep at a sufficiently large optical 
depth that the stellar light is fully absorbed, and the 
radiation diffusion approximation is valid. This choice of 
surface b oundary conditi on has also recently been advo- 
cated bv llro et, alJ l)2005l) . based on the results of time- 
dependent radiative models for the atmosphere of HD 
209458b. We motivate our choice with a simple toy prob- 
lem, and then discuss its relation to detailed radiative 
transfer solutions for the atmosphere. 

The atmosphere is heated by absorption of starlight, 
and possibly dissipation of day-night winds and tidal 
flows. Let there be an energy deposition rate e per unit 
volume in a radiative region of thickness Az. Choose 
boundary conditions T = at the top (for simplicity) 
and outward flux F = at the base of the heated layer. 
The latter choice is required in steady state so that the 
temperature deeper in the atmosphere not increase in 
time. The flux generated in the layer, which exits the 
planet, is F = sAz, and the temperature of the deep 
atmosphere is Td OC p ~ (tF/ct) 1 / 4 , where r = npAz is 
the optical depth, p is the density and n is the opac- 
ity. Hence an atmosphere subject to intense heating 
is expected to develop a deep isothermal region below 
the heated layer, the temperature determined primarily 
by the energy flux and depth of the layer, through r. 
This estimate of the deep isotherm temperature is simi- 
lar to that found for absorp tion of starlight for the proper 
choice of r l|Hubenv et al.11200^) . 

We now discuss the temperature profile for static atmo- 
spheres in more detail. In the absence of external irradi- 
ation, the photosphere of a planet will cool to a temper- 
ature T coo i ~ (-Fcooi/c) 1 / 4 ~ 100 K in a few Gyr's, where 
-Fcooi is the flux from the deep interior. A characteris- 
tic temperature at small optical depth for an irradiated 
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planet can be defined by balancing absorbed and emitted 
energy flux. For a star with mass M», radius i?* and ef- 
fective temperature T* a distance a = {GM^P^/Aii 2 ) 1 ^ 
away, this "equilibrium" temperature is 



o- — T* (_/?* /2o 



t 1400 K 



3 day 

P 



orb 



1/2 
1/3 , 



T, 



V6000 K 



/ 1 

Rq 



1/2 



Mq 

AL 



1/6 
(1) 



an order of magnitude larger than for an isolated planet. 
Hence the surface boundary condition is drastically al- 
tered from the isolated case. In general, the irradiated 
boundary condition will cause the planet to cool slower 
l)Burrows et al.ll2000l) . as we discuss in detail. As signifi- 
cant horizontal temperature variation is expected above 
the photosphere, T cq is an average temperature which 
gives the correct outgoing flux. Eq. assumes zero re- 
flection of the stellar photons, and should be multiplied 
by (1 — A) 1 / 4 for nonzero Bond albedo A. 

The (optical) incoming stellar photons not scattered 
back out of the planet are absorbed at the starlight's 
photosphere, typically at a pressure < 10 6 dyne cm~ 2 . 
Radiative balance implies an outgoing (infrared) flux 
F ~ (T C q/T cool ) 4 F cool ~ 10 4 F cool generated by ther- 
mal emission. This large flux may lead to a significant 
increase in temperature above T eq (the greenhouse ef- 
fect, e.g. Hubeny et al. 2003). This situation continues 
to a depth at which the starlight is fully absorbed, at 
which point the temperature profile becomes isothermal. 
Hence, a semi-infinite atmosphere subject to external ir- 
radiation, and with no internal flux deep in the atmo- 
sphere, becomes isothermal at large optical depth. We 
label the temperature of this deep isotherm Td ee p- Now 
including the internal cooling flux Fcoolj the temperature 
will again rise toward the interior, the gradient eventu- 
ally becoming large enough for convection to occur. 

Since the cooling luminosity is generated in deep layers 
with sufficiently large optical depth that the stellar light 
is fully absorbed, the radiation diffusion approximation 
is valid there. Furthermore, we will show in §0that the 
temperature profile becomes isothermal within a pres- 
sure scale height of the radiative-convective boundary. 
Hence the problem of determining the cooling luminos- 
ity is insensitive to many of the details of the absorp- 
tion of starlight. The only input needed from the full 
radiation transfer problem near the photosphere is the 
temperature of the deep isotherm, Td CC p- 2 

For tidally locked planets, the day side will be signif- 
icantly hotter than the night side in static atmospheres 
with negligible day-night winds. A more uniform temper- 
ature distribution results if winds can carry heat from the 
day to the night side without suffering radiative losses 
(e.g. Iro et al. 2005). We will show that the cool- 
ing luminosity is determined in deep layers with ther- 
mal time t t h ^ 10 3 yr. While significant day-night tem- 
perature asymmetries may exist near the optical photo- 
sphere, winds moving at even a tiny fraction, ~ 10~ 5 , of 
the sound speed could deposit heat on the night side in 

2 We expect that the degree to which this layer is isothermal 
depends on the number of pressure scale heights separating the op- 
tical photosphere from the radiative-convective boundary. Larger 
irradiation and lower core entropy should make this layer more 
nearly isothermal. 



less than a thermal time at the depths where the cool- 
ing luminosity is determined. Hence we have a strong 
expectation of a near-spherically symmetric, isothermal 
temperature profile deep in the radiative layer. 

3. NUMERICAL MODELS FOR THE INTERIOR 

In the deep interior where the diffusion approximation 
is valid, w e solve the mechanical a nd thermal structure 
equations (Landau & Lifshitz 1959) 



dm 2 
— — = Airr p, 
dr 

dP 

dr r 2 
dT _dPT^ 

dr dr P 



Gmp 



dl 
dr 



dm 
dr 



T 



ds 
~dt 



(2) 
(3) 
(4) 
(5) 



for the interior mass m, pressure P, temperature T, and 
outward luminosity I, as a function of radius r. Here S is 
the entropy per gram, and V = d In T/d In P is the loga- 
rithmic temperature gradient. The energy generation e is 
set to zero throughout this paper, as we study passively 
cooling planets. The subscript "cool" on the luminosity 
will be assumed for the rest of the paper. As the eddy 
turnover time is much shorter than the cooling time and 
convection is quite efficient, entropy is very nearly con- 
stant in space in the convection zone, but decreases in 
time due to cooling. Hence we treat dS/dt as a con- 
stant in the convection zone. For numerical convenience, 
we use this same value of dS/dt in the surface radiative 
zone. A negligible luminosity is generated there how- 
ever, so this error does not affect our results. While the 
entropy equation JSJ is valid on timescales longer than 
an eddy turnover time (~ yrs) in the convective core, 
the assumption of nearly spatially constant luminosity is 
only valid in the radiative envelope on timescales longer 
than the thermal time (~ 10 3 yr) there. As this is much 
shorter than the global cooling time, we expect our nu- 
merical cooling models to be as accurate as a relaxation 
(Henyey-type) code. 

The equation of state (EOS) from lSaumon et a 
(SCVH) is used with a mixture of 70% hydrogen and 
30% helium, ignoring metals, and using the tables which 
smooth over the plasma phase trans ition. There are been 
several improvements to SCVH (jSaumon et alJ 119991: 
iFortnev fc Hubbardll2004D using recent laser shock- com- 
pression data and including the effects of helium phase 
separation, which change the radii at the few percent 
level. We use the solar composit ion "condensed" phase 
opacities from lAllard et all l)2001t) . which includes the ef- 
fects of grains in the equation of state, but ignores their 
opacity, as is appropriate if the grains have condensed 
out. Mixing length theory is used to calculated V in 
convective regions, and radiative diffusion in radiative 
regions. The mixing length is set equal to the pressure 
scale height. 

Two boundary conditions are needed at the surface. 
First, the surface temperature is set to Tdeep, the temper- 
ature of the deep isotherm discussed in § [2j The second 
boundary condition is that we specify the surface to be 
at the (arbitrarily chosen) pressure P — 10 4 dyne cm~ 2 . 
As the surface layer is isothermal, the contribution to 
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the radius from near-surface layers is larger than for the 
radiative zero temperature profile, hence care is needed 
when comparing the radii computed here with previ- 
ous work. As the radius is somewhat dependent on 
the problem at hand (optical photosphere versus in- 
frared photosphere, corrections due to geometry in a 
transit, etc.) we have made this arbitrary choice of the 
surface for simplicity. The change in radius between 

pressures Pi and P 2 is AR — Jp 2 d\nP(kbT/fim p g) ~ 

(kbT / nm p g) ln(P2/Pi). For example, the radius must be 
decreased by AP = —0.022Pj for an outer boundary 
condition P = 10 6 dyne cm~ 2 for T = 1000 K and mean 
molecular weight fj, — 2.43 (70% molecular hydrogen and 
30% neutral helium). We do not include a solid core in 
the present calculations. 

We make a single model of a planet as follows. Planet 
mass M, core entropy S, and surface temperature Td eep 
are treated as fixed parameters. Assuming values for the 
planet's radius P, cooling luminosity L — l(R), dS/dt, 
and central pressure P c , we integrate outward from the 
center and inward from the surface. The four parameters 
are adjusted to make the integration variables (to, P, T, 
and I) continuous at a fitting radius. Given the subrou- 
tine to solve for a single model, evolving the planet in 
time is trivial. As we specify the core entropy S, and 
have solved for dS/dt, we compute the time it takes to 
cool from one entropy to the next. 

4. IMPACT OF IRRADIATION ON HEAT LOSS 

We now show the dependence of the cooling luminos- 
ity on the depth of the radiative-convective boundary, 
emphasizing the role of the opacity deep in the planet. 
Many of the luminosity dependences can be understood 
with purely local arguments, without the need to build a 
global planet model. Hence, qualitative statements can 
be made about cooling of EGP's under irradiation just 
given EOS and opacity tables. We make comparisons 
between the local arguments and the global numerical 
calculations as well. 

The convective core is capable of transporting enor- 
mous luminosities through fluid motion. Hence it is the 
large thermal resistance of the outer radiative envelope 
that determines the cooling flux. For an opacity which in- 
creases inward from the surface, this resistance is largest 
at the base of the radiative layer, hence it is the radiative- 
convective boundary that determines the cooling flux. 
This boundary is moved to higher pressures by irradia- 
tion ijGuillot et aUll99fl) . 

The outward flux carried by radiative diffusion is 

F = -^§, (6) 
inp dr 

where k is the Rosseland mean opacity. The maximum 
flux which can be carried by radiative diffusion is found 
using the adiabatic temperature gradient dT/dr | a( j = 
(V ad T/P)(-GTOp/r 2 ), where V ad = dlnT/dlnP| s (= 
2/7 for an ideal gas with five degrees of freedom) is the 
adiabatic temperature gradient. Multiplying by Airr 2 , 
the maximum luminosity per unit mass which can be 
carried by radiative diffusion at a local temperature T, 
pressure P, opacity k(T,P), and enclosed mass to ~ M 
is 

-V ad . (7) 
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Fig. 1. — Run of temperature vs. pressure for numerical models 
with T dcop = 500, 1000, 2000, 3000 K, Sm p /k b = 7, 8, 9, 10 and 
M/Mj = 0.32, 1.0, 3.2. Only the 19 curves in the age range 0.1 - 
10 Gyr are shown out of the total 48 curves. The circles show the 
position of the radiative-convective boundary. The triangles mark 
the center of the planet. Curves for different masses at the same 
S and Tdeep nearly overlie each other. The two nearly vertical 
lines show contours of constant density p = 0.01 g cm -3 (left) and 
0.1 g cm -3 (right). 



Choosing an entropy S, the right hand side of eq. J7J 
can be evaluated along an adiabat out from the center, 
yielding the cooling flux for a specified temperature T rc b 
at the radiative-convective boundary. Eq. J7J shows that 
the luminosity per unit mass depends solely on the en- 
tropy and irradiating flux, and that the luminosity is 
proportional to the planet's mass. 3 

Figure n snows the run of temperature versus pressure 
from numerical models for a range of M, S and Tdeep- 
Choosing S and T dee p, the temperature profile must fol- 
low the adiabat deep in the planet and the isotherm near 
the surface. An even stronger statement can be made, 
however. The temperature profile over the entire planet 
from the center to the top of the deep isotherm depends 
only on S and Tdeep, and is independent of M (since 
eq. {IT} and (O depend only on F/g oc L/M). Next, 
for a given irradiation flux (fixed Tieep)) the radiative- 
convective transition burrows deeper into the planet with 
time (decreasing S). Increasing the irradiation flux at 
fixed S also moves the radiative-convective region deeper 
into the planet. 

Temperature changes in response to small changes in 
flux in optically thick regions occur on the thermal time, 
estimated from eq. (j2J to be 



Uh = 



PC P T 



10 4 yr 



9F 

( 10 3 cm s~ 2 
V 9 



10 s dyne cm -2 
10 4 erg cm~ 2 s~ 



T 



10 3 K 



F 



(8) 



M 



kP 



3 It is commonly stated that the luminosity decreases with in- 
creasing mass. This is true at fixed core temperature, rather than 
fixed entropy. The derivatives can be related by d\nL/d In M\t = 
dXa.L/d\nM\ s + dlnL/dS\ u dS/dXa.M\ Tc . At fixed core temper- 
ature, entropy increases for decreasing mass (see Figure 151. 



5 



1Q29 =- 




1000 



2000 

T[K] 



3000 



Fig. 2. — Maximum radiative luminosity (scaled to M = lMj; 
see eq. Q) which can be carried by radiative diffusion vs. temper- 
ature along adiabats. The temperature T should be interpreted as 
T rc b, the temperature at the radiative-convective boundary. The 
five curves are for entropies Sm p /ki, = 6,..., 10, from bottom to 
top. 
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Fig. 3. — Luminosity vs. temperature of the deep isotherm from 
numerical calculation. Note that regions of positive slope in figure 
[5] correspond to zero slope here, implying the radiative-convective 
boundary jumps to a depth at which the slope is again negative. 
Calculation is for M = Mj, but can be extended to other masses 
using L oc M. The lines represent core entropies Snip/hf, = 6, 10 
from bottom to top. 



Here we have used typical numbers from Figure ^ for the 
radiative-convective boundary. Note that this estimate 
is much longer than the adjustment time near the optical 
photosphere (~ days, e.g. Iro et al. 2005), as the cooling 
flux is ~ 10 4 times smaller than the stellar flux, and the 
heat content increases oc TP. As the thermal time at 
the radiative-convective boundary is so much longer than 
the horizontal sound travel time (~ days), we expect the 
day-night temperature asymmetry to be small there. 
Figures and [3| show the local calculation of L/M 




S[k b /m p ] 



Fig. 4. — Luminosity per unit mass (scaled to 1 Mj) vs. en- 
tropy from numerical integrations. The lines correspond to sur- 
face isotherms T^eepfii"] = 500 (solid black), 1000, (dotted red), 
2500 (short dashed green), 3000 (long dashed blue), 3500 (dot 
short dash cyan). The lines for T(j ee p = 1500 and 2000 K over- 
lie those for T^ecp — 

1000 K (see Figure Three masses 

M/Mj = 0.32, 1.0, 3.2 are plotted, but closely overlie each other 
(except the largest masses at low entropies which are never reached) 
showing L/M is independent of mass at fixed S and T^ eep . 



evaluated along adiabats, and the global calculation of 
L/M versus Tdocp, respectively. The x-axis in Figure [21 
is the local temperature, which should be interpreted as 
T rc b, the temperature of the radiative-convective bound- 
ary. Care must be taken in Figure [2] in regions where 
L(T) increases inward. As we show in §0 L(T) must de- 
crease inward in order for convection to begin. Hence, if 
the chosen isotherm intersects a region of positive slope, 
such as the bump near T = 2000 — 2500 K, the convec- 
tion zone actually begins at a deeper point at which the 
slope L(T) is again negative. Such regions correspond 
to the flat parts of the curves in Figure [21 The result 
is that the luminosity generally decreases with irradia- 
tion temperature, or is roughly constant, but should not 
increase. This is t he origin of the resul t found by pre- 
vious investigators (Bu rrows et al]l20 00) that irradiated 
planets cool slower. Comparison of Figures [21 and show 
rough agreement in regions where L(T) is decreasing, the 
main discrepancies due to the ratio T rc b /Tdcop not being 
precisely a constant (see Figure QJ. 

Figure 0] shows luminosity versus core entropy for the 
numerical models. If Td eC p is constant during the evo- 
lution, Figure 0] shows the change in luminosity as the 
planet cools. Comparison of lines with different Td CC p 
clearly shows the monotonic decrease in luminosity as the 
irradiation temperature is increased. Aside from models 
with large mass (M — 3.2Mj) and irradiation tempera- 
ture (Tdocp = 3500 K) at entropies so low (S < 8kb/m p ) 
as to be unreachable in a Hubble time, the luminosity 
is proportional to the mass and the curves overlie each 
other. 

5. RADIATIVE-CONVECTIVE BOUNDARY 
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Fig. 5. — Solid lines show logarithmic temperature gradient V 
as a function of pressure for models with M = Mj, S = Sk^/nip 



and T, 



deep 



[A1 



500,1000,2000,3000 from left to right. The 



two dashed lines show the analytic formula in cq. 1131 with pa- 
rameters (Voo, a, P deep [dyne cm -2 ]) = (0.5,0.0,3.5 X 10 7 ) and 
(0.5,1.0,2.0 X 10 9 ). The upper envelope is set by the adiabatic 
gradient V a( j in the convection zone. Note the appearance of a 
radiative window which causes the radiative-convective boundary 
to be at P ~ 10 8 ' 5 dyne cm -2 for a large range in Td ee p- 



We now derive a simplified analytic model for the tem- 
perature profile at the transition from the surface radia- 
tive zone to the core convection zone. We relate Td CC p 
to T rc b, the radiative-convective boundary temperature 
where the cooling luminosity is determined. The scalings 
derived here are used in § to derive the scalings of the 
cooling luminosity. 

We assume constant gravity g, ideal gas pressure 

b — 



P = pkbT/fimp and power-law opacity 



KQ p a T 



K\P a T a . In the radiative zone, 

16erT 3 5 dT a + 1 16ag dT 4+a - b 



F = - 



3« dP 



b 3/ti dP^ 1 



(9) 



When integrating this equation, it's essential to retain 
the constant of integration. Defining the temperature 
gradient for a radiative zero solution V M = (a + l)/(4 + 
a — b) , we find 



T i+a - b = constant + V; 



3kiT 



P' 



a+l 



(10) 



At small pressure, T 
ture profile as 



16(7 g 

Tdeep, so we write the tempera- 



T = T dc 



1 + (P/P dccp ) a+1 



l/(4+o-6) 



which becomes isothermal below the pressure 



deep " 



wr c r b y /(Q+1) 

3 Kl F J 



(11) 



(12) 



4 Significant features in the opacity may be treated as broken 
power-laws. 



The logarithmic temperature gradient is then 

~i+(p/p deeP r +1 ' (3) 

which decreases sharply over a pressure scale height. A 
plot of V versus P is shown in Figure El for several val- 
ues of Tdeep- The upper envelope of the curves is set by 
the adiabatic gradient in the convection zone. Increas- 
ing Tdeep moves the boundary inward along the adiabat, 

Tdeep cx T d 1 ( ( C p ad , aside from regions where the opacity 
changes irregularly. Eq. (|13|) agrees well with the numer- 
ical integrations in regions where the opacity is smooth. 

To solve for the transition from radiative to convective 
zone, we set V = V a d- As a + 1 > 0, one must have the 
inequality > V a d for a convection zone to exist. We 
find the temperature and pressure at the boundary are 

^ l/(4+o-6) 

Tcb Tdeep 



Trcb — Tdc 



Voo — V a d 

Vad 
Vqo — V ad 



l/(a+l) 



so they differ by a factor of order unity from Tie 



(14) 
and 



Tdeep unless |V a d — Voo| << Voo- The decrease of 
V a d at large pressures seen in Figure makes the ra- 
tio Tcb/Ticcp closer to unity for large Tieep and T, as 
seen in Figure ^ Also note that at fixed S and Tieep, 
T rc b is largely independent of mass. 

6. ANALYTIC COOLING MODEL 

In § 01 we found that the luminosity scales with core 
entropy and irradiating flux over much of the relevant 
parameter space. Here we show that the core temper- 
ature also scales simply with mass and entropy when 
sufficiently degenerate. As a consequence, we derive an 
analytic model in which entropy has a simple power-law 
time dependence at late times. We compare the power- 
law model against numerical time integrations. 

The scaling of luminosity with Tieep and S can be 
found by substituting P(T,S) into eq. 0. Using the 
thermodynamic relation 

dS dT „ dP 



T 



'ad' 



P ' 



(15) 



and expanding about a reference point T re f , P re f and S le f , 
the adiabat is 

where AS = S — 5 rc f , and we approximate V a d and C' p 
as constants. For an ideal gas, C p V a d = kb/pm p , but 
particle interactions and molecular dissociation reduce 
CpVad below the ideal value. Inserting eq. ljlT)|. |Tl)l. 
and the power-law form of the opacity into eq. Q, we 
find 



T~T, 



Tieep \ a S rc {) 

exp (3— — 

T rc f / [ k b /m p 

where the exponents are (Figures |21 and 0} 

■o-6) f J^-l I ~ 0.0 -10.0, 

V Vad 

h/mp 



(17) 



a~(4 
/3~(a+ 1) 



Cp V a d 



2.5 - 3.5. 



(18) 
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Fig. 6. — Core temperature vs. mass. The five groups of lines 
show entropies Sm p /ki, = 6, 10 from bottom to top. Each group 
contains surface isotherms T dcep [K] = 500, 3500, inducing larger 
spread for small entropy. 



Fig. 7. — Core temperature vs. entropy. The three distinct lines 
are the masses M/Mj = 0.32,1.0,3.2 from bottom to top. For 
each mass, the irradiation temperatures T dccp [K] = 500,..., 3500 
are plotted, but there is so little dependence on T<j eep that the 
curves are indistinguishable. 



Examination of the exponent a shows that irradiation 
slows the cooling, i.e. luminosity decreases as Td ee p in- 
creases. The condition Voo > V a d is required for a core 
convection zone to exist, hence a > 0. Evaluation of 
a depends on the detailed density and temperature de- 
pendence of the opacity, which can be found in Figure |2 
Features of note are the positive slope near 2000 — 3000 K 
at which point a become small, and also the steep de- 
crease for T rcb > 3000 K. 

The exponent f3 can be estimated for an ideal gas (so- 
lar mixture, molecular hydrogen and neutral helium) and 
density independent opacity to be (5 ~ /i ~ 2.4. This 
ideal limit is expected for small Td CC p and hence low den- 
sity. As Tdccp is increased, molecular interactions make 
the gas less ideal, reducing the value of C p V a d and in- 
creasing p. This qualitative trend may be seen in Figure 

El 

The core temperature increases during the initial con- 
traction phase when the core is non-degenerate. A max- 
imum is reached when kbT c ~ Ep , the Fermi energy, and 
subsequently T c decreases as entropy decreases. In this 
degenerate phase, the core temperature depends mainly 
on mass and entropy, with only a weak dependence on ir- 
radiation. Figure shows core temperature versus mass 
for four adiabats and a range of irradiation temperatures. 
The dependence on Td ee p gives only a slight broadening 
of each adiabat. The dependence on mass is quite sim- 
ple when sufficiently degenerate. Figure \7\ shows the de- 
pendence of core temperature on entropy for a range of 
masses and irradiation temperatures, showing a simple 
exponential dependence at low entropy. For the degen- 
erate phase we write T c in the form (see Figure HjJ) 



T C (M, S) = T c<ie{ 



M 



of 



exp 



.(S — SVcf) 
h/mp 



(19) 



Using hydrostatic balance P oc M 2 /R 4 , and parameter- 
izing R oc M x , we estimate the exponents to be 

7 ~V ad (2-4A) -0.6-0.7 



S ~fcb/C p m p ~ 0.5. 



(20) 



Next we solve for the change in core entropy with time 
for the analytic model. We treat Td ee p and M as con- 
stants during the evolution. The entropy equation inte- 
grated over the convective core gives 



dS_ 
~dt 



L/M 

Jt7 7 



(21) 



where / = J(dm/M)(T/T c ) ~ 0.6 - 0.7 and we have 
treated dS/dt as constant in space. Plugging eq. Ijl9|l 
and 117|) into eq. H21fl . we find the following solution for 
the entropy with time 



exp 



S SVcf 



= I1 + T S 



-1/(0-5) 



h/m p 

where the characteristic cooling time is 
t s (M, T dcop ) 

/ \ ( kbT c /lTlp\ / Tdeep 



(22) 



/3-6 



L/M 



rof 



T, 



M 



(23) 



This solution has a number of notable features: 



(1) The fiducial evolution time, MkbT c /m p L, is 
just the initial time to radiate away the core's ther- 
mal energy. 

(2) At late times, t > tg, the solution is a power- 
law in time. As the initial cooling time is very 
short, the power-law occurs for most of the planet's 
lifetime. 

(3) The exponent of the power-law involves the 
change of luminosity and core temperature with 
respect to entropy. 

(4) Evolution timescale is slowed for large irradi- 
ation or large planet mass. The dependence on 
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Fig. 8. — Core entropy vs. age for a planet mass M/Mj = 0.32 
(dotted red), 1.0 (solid black), and 3.2 (dashed blue). For each 
mass, surface isotherms T^ ecp [K] = 1000,2000,3000 are shown 
from bottom to top. Luminosity is independent of 7d ee p for 
^deep = 1000 — 2000 K, hence the curves nearly overlie each other 
(see Figure 151. 
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Fig. 9. — Radius vs. mass curves for different core entropy and 
irradiation. Each group of lines with a different color and line 
style denote different entropies. Solid black, dotted red, and short- 
dashed green, long dashed blue and dot-dashed cyan lines represent 
entropies Snip/kt = 6, 10. Each group of lines represents deep 
isotherms T dcop = 500, 1500, 2500, 3500 K, from bottom to top. 



planet mass comes purely from the dependence of 
T c on planet mass (at fixed entropy). The depen- 
dence on irradiation is primarily through the lumi- 
nosity. 

Figure [S] shows entropy versus time for a range of mass 
and irradiation. Note the large spread in S at a given 
age. For low irradiation, S ~ 6 — 8kb/m p in the age 
range 1 — 10 Gyr, while S can be as high as ~ 10kb/m p 
for Tdecp = 3500 K. While the range of Td 00 p shown here 
gives fairly good power-laws, we note that at Td CC p < 
500 K there is a break occurring at ~ 1 Gyr, due to 
the increase in luminosity seen in Figure 0] below S = 
8kb/m p . 

While we have used the SCVH EOS and Allard 
et.al.(2001) opacities for numerical estimates, the ana- 
lytic solution makes it particularly clear which quantities 
need by evaluated for a given opacity table and EOS. The 
luminosity is sensitive only to the local conditions at the 
radiative-convective boundary, while the core tempera- 
ture involves building static (i.e. not time-dependent) 
models. 

7. RADIUS EVOLUTION 

Planets are initially nondegenerate in their core, and 
undergo rapid contraction until the core is degenerate. If 
T e ff is constant during the contraction, the energy equa- 
tion d/dt(-3GM 2 /7R) = -4irR 2 aT* s is solved to find 
the change in radius with time (see, e.g. Bildsten et al. 
1997) 



/?(./) ^/?, ( 



2/3 



/ 300 K 

V ^eff 



4/3 



1 Myr 



1/3 



(24) 



The core temperature T c ~ GMfxm p /Rkb oc i 1 / 3 is in- 
creasing during the non-degenerate phase, and reaches a 
maximum when k b T ~ Ep. For an ideal gas, k^T/Ep 
is a function only of entropy, so that the maximum tem- 
perature would occur at the same entropy for all planet 



masses. Coulomb interactions suppress the value of V a d 
below 2/5, so that if T c oc AT 4 / 3 and P c oc M 10 / 3 , the en- 
tropy at maximum temperature will increase a bit with 
mass. This can be seen in Figure as the two higher 
masses have maxima at higher entropy (off the plot) than 
the lowest mass. 

Once degeneracy sets in, the radius is primarily deter- 
mined in the degenerate core of the planet, although as 
irradiation is increased the contribution from the outer 
envelope becomes larger due to the increased scale height. 
This is clarified by writing the radius as an integral over 
pressure 



r(P) = I d\nP ( — 
lp \P9 



(25) 



In the degenerate core, P/p oc Ep while in the nonde- 
generate envelope P/p oc k b T. The contribution from 
the core is larger when Ep k b T unless the number of 
pressure scale heights in the envelope is much larger than 
the core. 

The equation of state in the core for M < Mj is com- 
plicated by strong Coulomb interactions. For illustrative 
purposes, an approximate equation of state including the 
leading order contributions from Coulomb interactions as 
well as ideal ion pressure is 

3 Z 2 e 2 \ pk b T 

(26) 



P 



10 cii 



where a, = (Anp/ piirip) 1 ^ 3 is the mean ion spacing and 
Piirip is the mean ion mass. The energy scales rele- 
vant for the core are Ep = (h 2 /2m e )(37r 2 p/ p e m v ) 2 l 3 ~ 
26 eV(p/i" 1 g cm" 3 ) 2 / 3 , k b T ~ 0.9 eV(T/10 4 K), and 



-Ecoui — Z e /oj 



20 eVZ 2 (ppig cm- 3 ) 1 / 3 is the 



Coulomb interaction energy between a nucleus and uni- 
form electron cloud. Ignoring the ion pressure term, the 
density at zero pressure is p zp ~ 0.2p e Z 2 g cm- 3 . As 



6 7 8 9 10 
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Fig. 10. — Radius vs. core entropy for different masses and 
irradiation. Solid red, dotted black, and dashed blue lines repre- 
sent masses M/Mj = 0.32, 1.0, 3.2. Each group of lines represents 
Tdoop = 500, 1500, 2500, 3500 K, from bottom to top. 
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Fig. 11. — Comparison of computed fractional radius (R(S) — 
Ro)/Ro (solid line) to fitting formula R(S) = Ro + 5R s e^ m P S ^ k b 
(dashed line) for a M = 1 Mj planet with T^eep — 

1000 K. 



the central density in Jupiter mass objects is near p zp , 
Coulomb interactions (and further the tendency to form 
bound states) stiffen the EOS, and are important in de- 
termining the radius. 

Figure [5] shows mass versus radius for a range of core 
entropy and irradiation. The effects of irradiation are 
seen to be most severe at low mass and low entropy, since 
Tdcop is becoming a significant fraction of the core tem- 
perature. At M ~ Mj/2 and low entropy, the range of 
irradiation temperatures shown here can change the ra- 
dius by as much as 50%. Radii for fully adiabatic planets 
(not shown here) agree well with the Td OC p = 500 K lines. 

Figure 1101 shows radius versus entropy for a range of 
masses and irradiation temperature. At late times in the 
evolution when the entropy is small, the radius is con- 
verging to some constant value which depends on both 
M and Td CC p- If the planet were allowed to cool under a 
constant irradiation field ind efinitely, it would approach 
an isothermal state (|Hubbardlll977|) at T = Td CC p with 
a radius 5 R = Rq. Although in practice planets will 
never reach this isothermal state, it is the minimum ra- 
dius to which the planet is evolving. Furthermore, it is 
the deviation around the isothermal radius, SR. = R— Rq 
which is changing with age. As we now show, SR has a 
particularly simple behavior with time over the entire 
observable range SR < R. 

To motivate the following numerical calculations, we 
first discuss the change in radius for a fluid element in 
mass shell m as the entropy is changed. The radius of a 
mass shell in the convection zone can be written 



o rr, 



dm' 



p(m', S) ' 



(27) 



5 In principle, Rq can be calculated by integrating the structure 
equations for a given EOS. In practice, such low temperatures and 
high densities are not covered by the SCVH EOS. In this paper, 
we compute the isothermal radius by fitting evolutionary curves of 
radius versus entropy, defining Ro by extrapolating to the small 
entropy limit. 



hence for fixed interior mass the change in radius with 
respect to entropy is 



dr 
dS 



1 



dm' dp{m',S) 



4irr 2 J p(m', S) dS 



(28) 



Given an equation of state p(P, S) , and switching to ra- 
dius as the integration variable, we find 



dr 1 / ,2 , , 

ds = -~ ' ' '"' 



1 ^p lp+T _ ld J^P 



CpdlnT 



dS 



where eq. ()15f) has been used. The second term in eq. 
H29|) mainly corresponds to a uniform shift in pressure in 
the core, due to the radius changing. Near the surface 
this term must go to zero since pressure is proportional 
to external mass, which is fixed. Consequently, the first 
term is most important. From eq. H26|l. the volume ex- 
pansion term is 01n/?/<91nT|p oc kbT/Ep, with a signif- 
icant correction due to Coulomb interactions which acts 
to increase the expansion since the electron pressure is ef- 
fectively lowered. Hence the change in radius in the core 
is proportional 6 to T c . As T c depends exponentially on 
the entropy (eq. |19|). the contribution to the radius from 
the degenerate core depends exponentially on entropy. In 
the nondegenerate envelope, din p/dlnT\p ~ —1. Plug- 
ging this result into eq. I|29|l implies that the change in 
radius due to the nondegenerate envelope scales linearly 
with entropy. As a consequence, it is less important than 
the exponential dependence from the core. 

A suite of evolutionary calculations has been done for 
M/Mj = 0.32, 1.0, 3.2 and T deep [K] = 500, 1000, 3500 
starting from high entropy and evolved to ages greater 
than 15 Gyr. Given the run of R(S), we fit a function 



R(S) = Rq + SRq exp(r)m p S / kb) 



(30) 



6 Eq. 1261 has ignored contributions to the Coulomb correction 
which depend on temperatu re, and do not scale linearly w ith tem- 
perature. Using the EOS in Potekhin & Chabricr (2000), we find 
the contribution of these terms to the volume expansion seems to 
be somewhat smaller than the ideal ion pressure. 
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TABLE 2 

Parameters for the Fi tti ng Function R(t) in eq. 



Fig. 12. — Fractional deviation in radius vs. core entropy for 
different masses and irradiation. This curves give a sense of how 
much the radius changes during evolution. Solid red, dotted black, 
and dashed blue lines represent masses M/Mj = 0.5, 1.0, 1.5. Each 
group of lines represents T^cep = 500, 1000, 3500 K, from left to 
right. Only ages in the range 0.1 to 10 Gyr are shown for each 
curve. 



to determine the isothermal radius Ro, coefficient SRo, 
and exponent 77. The coefficients Rq, SR and r\ depend 
on M and T^ ccp . The small entropy points were more 
heavily weighted to force the fit to agree there. The 
weighting was adjusted until the fit agreed for as large a 
region in S as possible (for the plots here we used weight- 
ing oc R 10 .) A comparison of the fit against the data for 
one example is given in Figure 1111 The agreement is 
good at small entropies, and gets worse for large entropy 
as degeneracy is lifted. We find good agreement between 
7/ ~ 0.5 — 0.7 and S from eq. (jT§|) . as expected if 6R oc T c . 

The deviation of the radius about the isothermal value 
is plotted for all runs over the age range 0.1 — 10 Gyr 
in Figure ^| Recall that Ro is different for each line. 
Note that each line is approximately a power-law, even to 
SR/Rq ~ 1, where the degenerate approximation breaks 
down. Hence the fitting formula often works better than 
naively expected. 

We now combine the power-law cooling result in eq. 
(J22J) and with the fit for the radius in eq. ||2U|) to 
find 

-n/{0-S) 

(3.1) 



t 



R(t) = R + SR eyip(r)m p S Te {/k b ) 1 + 

V *S 

At late times t 3> t$, the deviation in radius from the 
isothermal value is a power-law in time. In order to pro- 
vide useful fits to our evolutionary tracks, we parametrize 
this late time power-law as 



R(t) = R a + 8R X 



1 Gyr 



r,/(0-S) 



(32) 



where Ro is again the isothermal radius and SRi = 
SRoexp(rjm p S le {/kb)(ts/l Gyr) ?) /^~' 5 ) is the deviation 
at an age of 1 Gyr. We fit tracks of R(t) to find the 
coefficients Ro, 8R1 and Tjjip — S) in the same way as 
the fits for R(S) in eq. (|30(l . The coefficients are given 
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Fig. 13.— Radius vs. age for M/Mj = 0.32, 1.0, 3.2 and T deep = 
500, 1500, 2500, 3500 K from bottom to top within each panel. Solid 
lines show the numer ical evolutions, while dashed lines show the 
fitting formula in eq. 1321 . 



in Table [5] Comparison between the numerical evolu- 
tionary tracks for R(t) and the analytic fit in eq. I|32[l 
are given in Figure ITTfl The agreement is generally very 
good. 

Approximate values and scalings of the coefficients in 
Table [5] can be understood as follows. The expected 
power-law index r//(f3 — S) ~ 0.6/3.0 = 0.20 agrees well 
with the temperature range Td CC p = 1000 — 2000 K where 
L is independent of Td ccp . At large irradiation, Fig- 
ure El shows r\ increases and Figure 0] shows that (3 
decreases, explaining the increase in rj/{[3 — 5). Since 

SRi cx tg^ -15 ' oc ~ , regions of constant (de- 

creasing) slope in Figure [21 correspond to SR\ being con- 
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Fig. 14. — Td CC p vs. age for planets with a given M and i?. 
Shown are lines of constant radius for a given mass, labeled by 
(M/Mj, R/Rj). At a given age, this plot shows the value of T deep 
required to explain a certain mass and radius. There is significant 
degeneracy between Td ecp and age. 



stant (increasing). The magnitude of 8R\ can be esti- 
mated from FigurelT^land eq. (|23|) . Interestingly, Rq can 
be somewhat bigger for M = 0.32A// than for the higher 
masses. While a larger radius is expected for strong irra- 
diation, we caution the reader about interpretation of the 
exact values for Rq. It would be interesting to compare 
the values obtained by fitting tracks with actual calcula- 
tions of isothermal planets given a sufficiently accurate 
low temperature EOS. 

Given measurements of planetary mass, radius and 
age, Tdeep can be constrained. Figure PHI shows the value 
of Tdeep required to explain a planet of a given mass and 
radius, as represented by different lines, as a function of 
age. The lines slope up to the right since the cooling 
must be slower (higher Tieep) to reach the same radius 
at larger age. As the lines are not horizontal or vertical, 
there is significant degeneracy between Tieep and age. 
Large radii in the age range 1 — 10 Gyr can only be ex- 
plained by large irradiation temperatures for the mass 
range 0.32 — 3.2Afj. For each mass and radius, there is 
a minimum age which is set by the unirradiated planet, 
resulting in a steep slope down to the left. We shall use 
Figure ITU in § [Slto constrain Tj ee p for the observed tran- 
siting planets. 

8. APPLICATIONS TO TRANSITING PLANETS 

We now compare our theory to the observed masses 
and radii of the transiting planets (Table 0). Figure 
1151 shows radius versus mass for the observed tran- 
siting planets. The points with errorbars are the 
data. The three different hatched regions show Td OC p = 
500, 2500, 3000 K from bottom to top. The change is 
gradual from Tieep = 500 to 2500, and then acceler- 
ates for higher temperatures (see Figure[2J). Within each 
hatched region, a spread of ages from 1 (top) to 10 Gyr 
(bottom) is shown. The radius of HD 149026 is so small 
as to be well outside the plot. It clearly has a large 
abundance of heavy elements. The radii of the other 
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Fig. 15. — Radius vs. mass for observed transiting planets 
compared to our cooling model. The points with error bars give the 
observed masses and radii, as listed in Table The lower (red), 
middle (green), and upper (blue) hatched areas denote T^ eep = 
500, 2500, 3000 K. For each T dee p and M, the spread in R denotes 
ages 1 Gyr (top) to 10 Gyr (bottom) in each hatched region. HD 
149026's small radius places it well outside the plot. 



eight planets can be broadly explained with solar com- 
position, ages in the range 1 — 10 Gyr, and temperatures 
deep in the atmosphere Td CC p < 3000 K. The largest 
radii requiring the most irradiation are HD 209458, HD 
189733 and OGLE-TR-56. 

There are significant uncertainties in fitting stars on 
the main sequence to find stellar ages. Hence there is 
motivation to understand how a range of ages affects 
the range of observed radii. Figure ^| shows devia- 
tion from the isothermal radius by factors 1.1 — 2 in 
the age range 0.1 — 10 Gyr. The length of each track 
gives an idea of the uncertainty in radius due to an 
uncertainty in age. Using the fitting formula in eq. 
<|32|) . the fractional difference in radius between ages 
ti and t 2 is ~ [7?/(/3 - 8)]{5Ri/R Q ) hx{t 2 /tx). For the 
strongly irradiated case, if we choose characteristic val- 
ues r//(f3 — 5) = 0.5, SRi/R — 0.5, and a factor of two 
error in age ti — 2t\, the fractional error in radius is 
17%. Hence, the age dependence for strongly irradiated 
planets is important because (i) the decrease in time is 
steep, and (ii) strong irradiation increases the size of <5i?i 
relative to <5i? - 

Next, the parameter Tieep is crucial for the cooling 
rate, but is not directly measurable. Here we constrain 
Tieep using measured mass, radius and age. We then 
compare Tieep to the equilibrium temperature. 

We interpolate over the age — Ti CC p tracks in Figure 
1141 for the mass and radii appropriate for each planet 
(except HD 149026, which we do not discuss). Since the 
uncertainty in Tieep due to the error bar in planet mass 
is smaller than that due to the error bar in radius, we fix 
the planet mass at the central value and only vary the 
radius. For those planets with an age range quoted in the 
literature, we show the age range in the plot, and derive 
the range of Tieep consistent with the age range. These 
values are listed in Table ^ F° r those planets with no 
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Fig. 16. 



Constraints on T dee p for HD 209458 and TrES- 



1. Lines labeled by (M/M.j,R/Rj) given allowed values of Tdecp 
versus age for the central value of mass, and a range of radii given 
by the radius error bar. Vertical lines give age range from main 
sequence fitting of parent star. Shaded region shows range of T dccp 
allowed given uncertainties in mass, radius and age. 
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Fig. 17.— Constraints on T deep for OGLE-TR-56 and OGLE- 
TR-10. See Figure^llfor descript ion. 



age determination, we find the maximum value of T^ eep 
consistent with an age less than 10 Gyr. 

The constraints on Td ee p for all planets except HD 
149026 are shown in Figures dEI El and d and sum- 
marized in Tabled The T dc op of HD 209458 is best con- 
strained due to the small error bar on mass and radius, 
as well as detailed fitting of the parent star to find the 
age. From FigureQlwe find T dccp = 2200 - 2800 K; HD 
209458b is not consistent with an un- irradiated planet. 
OGLE-TR-56 has a weak lower limit on Tdccp which is far 
less than the equilibrium temperature. All other planets 
with age constraints have only upper limits, set by the 
upper limit on the radius, since the lower limit on the ra- 
dius is consistent with no irradiation. Plots are provided 
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Fig. 18.— Constraints on T dccp for OGLE-TR-111 and OGLE- 
TR-113. See Figure [TBI for description. 
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Fig. 19.— Constraints on T dccp for OGLE-TR-132 and HD 
189733. See Figure HSl for description. 



for those planets with no age constraints at the present 
time. 

The equilibrium temperature T oq = T ae (R*/2a) 1 / 2 is 
measurable, but plays no part in our model. On the 
other hand, the temperature of the deep isotherm is not 
measurable, but is crucial for the cooling rate. From 
radiative transfer models, we expect these two tempera- 
tures to be roughly proportional, the exact ratio deter- 
mined by the size of the greenhouse effect (Q- Hence, 
they should be strongly correlated. Figure 150] shows a 
plot of T eq versus Tdccp- The large error bars on Td ee p, 
due to large error bars on the radius, prevent one from 
drawing robust conclusions. It is in principle possible for 
a correlation (sloping up to the right) to exist given the 
current error bars, however, it is not required. 

Tighter constraints on Td ee p require the following: 

• Significantly smaller error bars on the radii of 
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• These scalings allow us to derive an analytic model 
for the cooling, which shows power-law decrease 
over a large range of parameter space. The part 
of the radius which changes in time (the deviation 
from the isothermal planet) is also a power-law. An 
analytic formula for radius evolution is given in eq. 
(|3"2l . with coefficients in Tabled 



• While we have used the SCVH EOS and Al- 
lard et.al.(2001) opacities for numerical estimates, 
the analytic solution makes it particularly clear 
which quantities need to be evaluated for a given 
opacity table and EOS. The luminosity is sensi- 
tive only to the local conditions at the radiative- 
convective boundary, while the core temperature 
involves building static (i.e. not time-dependent) 
models. 



Fig. 20. — Equilibrium temperature vs. Td eep . The temperature 
of the deep isotherm T(j eep is expected to be strongly correlated 
with the equilibrium temperature T eq , giving a line sloping up to 
the right. Bracketed lines give constraints on T^ ccp , while arrows 
pointing downward represent upper limits found when no age was 
available in the literature. 



We have compared our theory to observed masses and 
radii for the transiting planets in Tabled (except for HD 
149026, which clearly has a large abundance of heavy 
elements). Our findings are as follows: 



OGLE-TR-56 and OGLE-TR-132. 

• An age estimate is needed for HD 189733. If it 
is found to have an age > 1 Gyr, Td ee p will be 
well constrained with a value much larger than T eq , 
similar to HD 209458b. 

• Age estimates are needed for OGLE-TR-10, 
OGLE-TR-111, and OGLE-TR-113. However, 
given the present error bar on radii of OGLE-TR-10 
and OGLE-TR-113, Td eC p will be constrained only 
at the factor of two level. OGLE-TR-111 is an in- 
teresting case, as it must be older than ~ 3 — 4 Gyr 
to be consistent with our model. 

We encourage efforts in these directions. 

9. CONCLUSIONS 

We have presented calculations of cooling and radius 
evolution for strongly irradiated planets. Novel aspects 
of this model are the following: 

• We argue that the generic outcome of strong sur- 
face heating, whether it be due to absorption of 
stellar flux or dissipation of winds and tidal flow, 
is that a deep isothermal region exists above the 
radiative-convective boundary. The thermal time 
in this layer is sufficiently long that the tempera- 
ture profile is approximately spherically symmetric, 
irrespective of the size of the asymmetry near the 
photosphere. We assign this region the tempera- 
ture Tdccp and treat it as a boundary condition for 
the cooling models. 

• We show that the cooling flux is determined at 
the radiative-convective boundary, which is much 
deeper than the photosphere. Scalings of the flux 
with core entropy, Td ee p, and mass are computed. 



• Figure El shows mass versus radius for eight tran- 
siting planets, compared to our model. The radii 
can be broadly explained with solar composition, 
ages in the range 1 — 10 Gyr, and temperatures deep 
in the atmosphere Td eC p < 3000 K. The largest 
radii requiring the most irradiation to explain are 
HD 209458, HD 189733 and OGLE-TR-56. 



• Figures El El El an d El show constraints on 
Tdeep using measured masses, radii, and ages (when 
available), and their uncertainties. We find that 
only HD 209458b is well constrained, with Td 0C p = 
2200 - 2800 K. OGLE-TR-56 has a weak lower 
limit, and the other six planets have only upper 
limits, due to the large measurement uncertainty in 
the radius, or lack of an age determination. These 
constraints are summarized in Tabled 



• The equilibrium temperature T cq is measurable, 
but plays no part in our model. The deep isother- 
mal temperature Td 0C p is not measurable, but is 
crucial for the cooling rate. Radiative transfer cal- 
culations find these two temperatures should be 
strongly correlated. Figure El shows T cq versus 
Tdeep- As only upper limits on Td CC p are available 
for all but HD 209458b and OGLE-TR-56, it is dif- 
ficult to draw conclusions at the present time. It is 
in principle possible for a correlation to exist given 
the current error bars, however, it is not required. 

We hope that our models have illuminated the need for 
more accurate ages and radii. Once those are in hand, 
our calculations will provide a measurement of Td ee p of 
adequate accuracy to compare to T eq , thus constraining 
greenhouse physics and day-night transport. 
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